NearestNeighbor Subroutine

SUBROUTINE NearestNeighbor &
!
(network, grid, weights, gridPolygons)

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid

!Arguments with intent (out):
REAL, OPTIONAL, INTENT(OUT), POINTER      :: weights(:)
TYPE (grid_integer), OPTIONAL, INTENT(OUT) :: gridPolygons

!Local declarations
TYPE (ObservationalNetwork) :: activeNetwork
INTEGER (KIND = short) :: count
INTEGER(KIND = short), ALLOCATABLE  :: polygons(:,:)
REAL (KIND = float) :: dist, dist_min
INTEGER (KIND = short) :: i, j, k

!----------------end of declarations-------------------------------------------

!check for available measurements
CALL ActualObservations (network, count, activeNetwork)

!set points
point1 % system = grid % grid_mapping
point2 % system = grid % grid_mapping

!allocate polygons matrix
ALLOCATE (polygons (grid % idim, grid % jdim) )
polygons = -1


DO i = 1, grid % idim
	col:DO j = 1, grid % jdim
		IF (grid % mat(i,j) == grid % nodata) THEN
			CYCLE COL
		ELSE
			!initialize
			grid % mat(i,j) = 0.0
		END IF
		dist = -9999.99
		dist_min = -9999.99
		
		!compute distance from cell center to observations
		CALL GetXY (i,j,grid,point1 % easting,point1 % northing)
        DO k = 1, activeNetwork % countObs 
           point2 % northing = activeNetwork % obs (k) % xyz % northing  
           point2 % easting = activeNetwork % obs (k) % xyz % easting
           dist = Distance(point1,point2)
		   !search for nearest points
			IF ( dist_min == -9999.99 .OR. dist_min > dist ) THEN
					polygons(i,j) = k
				    dist_min = dist
			END IF
		END DO
	END DO col
END DO

!finalize interpolation
DO i = 1, grid % idim
	DO j = 1, grid % jdim
		IF (grid % mat(i,j) == grid % nodata) CYCLE
		k = polygons (i,j)
		grid % mat(i,j) = activeNetwork % obs (k) % value
	END DO
END DO

IF (PRESENT (gridPolygons)) THEN
    CALL NewGrid (gridPolygons,grid)
	gridPolygons % mat = polygons
END IF

IF (PRESENT(weights)) THEN
    ALLOCATE(weights(activeNetwork % countObs))
    weights = 0.
	    DO i = 1, grid % idim
		    DO j = 1, grid % jdim
			    IF (grid % mat(i,j) == grid % nodata) CYCLE
			    weights(polygons(i,j)) = weights(polygons(i,j)) + 1
		    END DO
	    END DO
    weights = weights / CountCells(grid)
END IF

RETURN
END SUBROUTINE NearestNeighbor